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We use the two-step density-matrix renormalization group method to elucidate the long-standing 
issue of the universality class of the Mott transition in the Hubbard model in two dimensions. We 
studied a spatially anisotropic two-dimensional Hubbard model with a non-perfectly nested Fermi 
surface at half-filling. We find that unlike the pure one-dimensional case where there is no metallic 
phase, the quasi one-dimensional model displays a genuine metal-insulator transition at a finite 
value of the interaction. The critical exponent of the correlation length is found tohe u ~ 1.0. This 
implies that the fermionic Mott transition, belongs to the universality class of the 2D Ising model. 
The Mott insulator is the 'ordered' phase whose order parameter is given by the density of singly 
occupied sites minus that of holes and doubly occupied sites. 



I. INTRODUCTION 

In the studies of the Mott transitioiiiiS, in the ground 
state of the Hubbard model^, there are well controlled re- 
sults in the pure one-dimensional (ID) case^ and in the 
limit of infinite dimensions^^^ only. In ID, there is no 
metallic phase, the Mott gap opens as soon as the in- 
teraction U > 0. In infinite dimensions, the dynamical 
mean-field theory which is exact predicts a Mott tran- 
sition at the critical coupling, Uc ~ W, W is the band 
width. However, the transition has mean-field critical ex- 
ponents. This anomaly is due to the local nature of the 
infinite dimensional solution. Hence, the one-dimensional 
and the infinite dimensional solutions may not be directly 
applicable to experiments. Studies of the Mott transi- 
tion in the Hubbard beyond these special limits of one 
dimension and infinite dimension are thus of crucial im- 
portance. 

For more than a decade, a great deal of effort has been 
devoted to applying quantum cluster theoriesSrii to the 
study of the Mott transition in the Hubbard model in 
two dimensions (2D). Quantum cluster theories include 
non-local correlations. They predict a finite critical value 
for the interaction at the transition. This critical value 
depends on the cluster size. However, when applied to a 
finite dimensional model, they are exact only in the limit 
of infinite cluster size. In quantum cluster theories, the 
effect of the interaction on physical quantities such as the 
single-particle Green's function is restricted to the clus- 
ter sites. The correlation are fully accounted for distances 
which are smaller than the cluster length, r ^ Lc- When 
r ^ Lc, the Green's function has an effective mean-field 
decay. Restricting the effect of the interaction at dis- 
tances r ^ Lc is probably justified away enough from 
the critical point where the correlations are expected to 
be short-ranged. A consequence of this restriction of the 
correlations to the cluster length is that the exponents 
at the transition are always mean-field like for a fixed 
cluster sizelS. A systematic finite cluster size analysis is 
therefore necessary for a correct description of the tran- 
sition. However, most of applications of quantum cluster 
simulations have been done on relatively small clusters. 



These are not enough to reliably predict the low-energy 
physics at the quantum critical point. 

Unlike the fermionic model, in the 2D Bose-Hubbard 
model which displays a transition from a superfluid to 
a Mott insulator, analytical approache a^^d^ and large 
scale Monte Carlo simulationsii have yielded reliable in- 
formation about its critical behavior. The transition for 
fixed boson density belongs to the universality class of 
the classical three-dimensional (3D) XY model. This has 
also been reported on the 2D Jaynes-Cummings-Hubbard 
modelA^. Unfortunately, for the fermionic Hubbard 
model Monte Carlo simulations predict Uc = 0. This 
is because of the nesting induced Slater transition^^'^''. 
In absence of perfect nesting, the Monte Carlo method is 
hampered by the sign problem. Large scale simulations 
are not possible. 

Recent interest has been raised by slave rotor 
analyses^ii^Si. These analyses suggest that the transition 
in the 2D fermionic Hubbard model may belong to the 
3D XY universality class as the bosonic Hubbard model. 
In Ref j M^' — a slave rotor representation of the fermionic 
operator Cia = hfia-, where bi is a spinless boson and 
fia- a charge-less spin, was used to map the Hubbard 
model to a free spinon Hamiltonian self-consistently cou- 
pled to a bosonic term (or XY term in a spin representa- 
tion of bosons) . The fermionic Mott transition is in this 
form a transition between condensed (Fermi liquid) and 
non-condensed (Mott insulator) phases of bosons. This 
factorization may be justified in the Mott phase where, 
because of the Mott gap, spin and charge degrees of free- 
dom may be separated. However, as the critical point 
is approached, is the gauge field weak enough to justify 
the decoupling between spin and charge? If not would 
that modify the critical behavior predicted by the slave- 
rotor approximation? Only a non-biased calculation of 
the Hubbard model can yield the answer. 

The slave-rotor prediction is in disagreement with an 
earlier approximate mapping^"^ of the Hubbard model 
to a generalized Blume-Emery-Griffiths modeP'* of the 
— mixtures with an additional term whose effect 
on the nature of the transition is not known. In this 
mapping, doubly occupied and empty sites corresponds 
to sites and singly occupied sites to sites. This 
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mapping suggests instead that the Hubbard model is in 
the universahty class of the Ising model. But the ex- 
tra term which accompanies the Blume-Emery-GrifSths 
model could well lead to another universality class. 

In a recent paper—, we reported a two-step density- 
matrix renormalization group (DMRG)^ study of the 
Mott transition in the ground state of the quasi-one- 
dimensional (ID) Hubbard model at half-filling. We find 
that in contrast to the pure ID case for which there is 
no metallic phase, there is an authentic Mott transition 
in the quasi-lD model. However, it is possible to argue 
that in the quasi-lD dimensional Hubbard model stud- 
ied in Refi^, the Fermi surface is perfectly nested, thus 
our analysis which predicts a gapless phase in the weak- 
coupling regime, would miss an exponentially small gap 
A oc exp — that would open as a consequence of a 
Slater transition. However, our numerical data did not 
support the existence of such a gap. Arguments support- 
ing a gap opening induced by perfect nesting are pertur- 
bative: the divergence of the non-interacting susceptibil- 
ity xo (q) at the nesting wave vector leads to that of the 
interacting spin susceptibility, Xs{<l) « 1/(1 ^ C^Xo(q))- 
However, the actual susceptibilities and interaction in 
the expression of Xs(q) are renormalized. Attempts to 
compute the renormalized susceptibilities and interaction 
within the self-consistent parquet formalismSI lead to in- 
tractable equations. Hence, the effect of these renormal- 
ization effects on the mean-field solution remains an open 
problem. 

In this paper, we present a well controlled study 
of the Mott transition in the Hubbard model with a 
non-perfectly nested Fermi surface beyond the special 
cases of ID and infinite dimensions. The choice of 
the non-perfectly nested Fermi surface precludes the 
theoretical possibility of a gap induced by the Slater 
anti-ferromagnetism mechanism. The two-step DMRG 
method is first checked on the transition between a para- 
magnetic and an anti-ferromagnetic ground states in the 
quasi-lD Heisenberg model with S = I. In agreement 
with a quantum Monte carlo study^^, we find that this 
transition belongs to universality class of the 3D classical 
Heisenberg model. For the quasi-lD Hubbard model, we 
find that, in contrast to the pure one-dimensional model, 
there is a genuine ground-state Mott transition at a finite 
critical value of the interaction. Data analysis of the crit- 
ical behavior of this model show that, in agreement with 
the mapping to the Blume-Emery-Griffiths model^^, the 
Mott transition in the 2D Hubbard model belongs to the 
universality class of the 2D Ising model. 



II. MODEL 




FIG. 1: The anisotropic frustrated lattice with longitudinal 
tj;, transverse ty, and diagonal td hopping parameters. 
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the hopping parameters tx, ty, and td, respectively in the 
longitudinal, transverse, and diagonal directions, are il- 
lustrated in FigHJ The presence of td ensures that the 
non-interacting Fermi surface is not perfectly nested, ty 
and td must be {ty , td) <C t^ for the two-step DMRG 
method to be accurate. In this study, we set t^ = I and 
ty = td = Q.Ohtx- The choice of this model thus pre- 
cludes the theoretical possibility of the nesting induced 
exponentially small gap. The band- width is W — '^Atx, 
we set u — U/W. 



III. TWO-STEP DENSITY-MATRIX 
RENORMALIZATION GROUP 

The two-step DMRG is a generalization of the conven- 
tional DMRG method^S to quasi-lD Hamiltonians. The 
DMRG is a RG procedure in which the reduced density- 
matrix is used to retain the most important states of 
the system. The DMRG itself is a crucial improvement 
over the block RG method^ which extended the Wilson 
RG method^i used in the solution of the Kondo impu- 
rity problem to lattice models. The block method has a 
major handicap, by dividing the lattice into independent 
blocks, it neglects at its initial step the inter-block inter- 
action. But if the inter-block interaction is of the same 
order as the intra-block interaction, this introduces an er- 
ror from which it is difficult to recover even by keeping a 
large number of states. In the DMRG the lattice is built 
by initially coupling the block to the rest of the lattice. 
Let us consider a system (S) coupled to an environment 
(E), let Ng and be respectively the number of states 
respectively of the system and for the environment. Let 

be for instance the ground-state wave function of the 
super-system including the system and the environment. 



We consider the Hubbard model with the local inter- 
action U and the following non-interacting single-particle 
energies. 



€{kx, ky) 



2tx COS kx 



2ity COS ky 



^{S,E) 



(2) 



where the ^/'i^'s represent the system's basis states and 
the 's the environment basis states; Ns and are re- 
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spectively the total number of states of the system and of 
the environment. The essence of the RG procedure is the 
truncation of the Hilbert's space, starting with a smah 
system for which the total number of states can be kept, 
at some step when the lattice gets large, only a smaller 
number nis < Ng of the system's states can be kept. The 
error in this truncation is given by the eigenvalues A^^ of 
the reduced density-matrix of the system, 

Ds= HS,EW{S,E). (3) 

From the relation, 

E ^'==1' (4) 

the error made by representing the system by mg states 
instead of Ng is given by, 

p = ^- E (5) 

For a large number of ID models, p is very small if rus is 
only a few hundreds. Application of the DMRG method 
to Heisenberg chains with 5* = 1/2 or 5 = 1—, rus ^ 
100, the ground-state energy, correlation functions and 
lowest excitation gap were obtained with an astonishing 
accuracy. 

It was hoped that, given the level of accuracy of the 
DMRG for ID models, the method would also perform 
reasonably well for 2D models. However, for a 2D lat- 
tice, the value of necessary to retain good accuracy 
appears to increase exponentially with the system size. 
This is related to the entropy area law which predicts an 
exponential increase of rris oc 2^ in 2D. The entropy 
area law implies that the direct application of the 2D 
DMRG would only be limited to relatively narrow sys- 
tems, it however leaves a window of success for quasi- ID 
systems as we will explain below. The study of quasi- 
ID models would yield valuable information about the 
corresponding isotropic models. Most importantly the 
two-step approach had a direct relevance to the physical 
properties of quasi-lD materials for which ty <^ such 
as the organic and inorganic quasi- ID conductors. 

Let us consider for instance the Hubbard chain with a 
charge gap A. If the transverse coupling ty is infinitely 
small with respect to A, so that the system remains in 
the same phase as the decoupled chains. It is obvious 
that the decoupled chain limit is a good starting point 
to describe the weakly-coupled chain system. As ty in- 
creases, the quality of decoupled chain as a starting point 
will decrease, if the same number of states is kept, un- 
til ty reaches a quantum critical point ty at which the 
systems enters in the 2D regime. In principle, when ty 
is in the 2D phase, it would be wrong to start from the 
decoupled chain limit. This is because there are a huge 



number of low-lying states with nearly equal weight in 
the reduced density-matrix. 

The important point which nevertheless makes calcula- 
tions possible is that actual calculations are done on finite 
systems which have a discrete spectrum. Thus even if ty 
has a value corresponding to the 2D phase for a system 
size L, given the discreteness of the energy spectrum for 
a finite system, if the energy width of the states kept 
is such that AE ty, starting from decoupled chain 
might still lead to accurate results. For such a system, 
the DMRG can be used to study the ground-state phase 
transition since it will display a different scaling behav- 
ior above and below ty. The same type of analysis may 
be used for gap-less chains as well, A(L) will yield the 
relevant energy scale above and below the transition. 

The separation of the energy scales is basic idea of the 
two-step DMRG^. The two-step DMRG uses the ex- 
traordinary accuracy that the DMRG can achieve in ID 
in two steps. In the first step, the low-energy Hamilto- 
nian is obtained accurately using the DMRG. Then, in 
the next step small transverse perturbations are inserted. 
The 2D effective Hamiltonian is ID, the DMRG is again 
applied to solve the problem in the transverse direction. 
Indeed, this procedure is valid only if the transverse cou- 
plings are very small with respect to the longitudinal cou- 
plings. The success of the two-step DMRG in yielding 
reliable results on the eventual new physics induced by 
the perturbation will depend on the value of the critical 
transverse coupling necessary to drive the systems in a 
new phase. If the magnitude of the perturbation ty nec- 
essary to drive the system away from the ID physics is 
small in comparison with the width of the states kept, 
the two-step DMRG is expected to be successful. This is 
for instance the case of coupled Haldane chains studied 
in section ITVl However, if the magnitude of the perturba- 
tion is too large, the two-step DMRG would not be able 
to describe the 2D physics accurately. 

The real challenge in the two-step starts after finishing 
making the program code work. The essential part of the 
subsequent activity is finding a region in the parameter 
space of a given model where interesting physical results 
can be extracted. For more details about the two-step 
DMRG, we refer the reader to ReiM. 

In the first step of the DMRG, we targeted charge sec- 
tors with Ne, Ne±l, Ne±2, where corresponds to the 
number of electrons at half-filling; for each charge sector, 
we targeted the spin sectors with the lowest Sz,Sz ± 1; 
hence we targeted a total of ritarg = 17 charge-spin sec- 
tors during each DMRG iteration. The reduced density- 
matrix was given by, 

Ds= E E '^k{S,E)<t>US,E). (6) 

k—l.ntavg ie—l-Ne 

where we assigned an equal weight Uk — 1/17 to each 
state In all the simulations we kept msi = 512 states 
such that the largest truncation error was pi ~ 10~^ for 
systems of up to = 32 as can be seen in Table HI 
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12 X 13 


16 X 17 


20 X 21 


24 X 25 


28 X 29 


32 X 33 


AE{u = 0) 


1.6220 


1.2683 


1.0410 


0.8819 


0.7685 


0.6772 


Pi(m = 0) 


8 X lO"'-' 


3 X 10" ' 


7 X 10"' 


1 X 10"" 


3 X 10"" 


4 X 10"" 


P2{u = 0) 

















4 X 10"* 


AE{u = 0.4261) 


1.5630 


1.2333 


1.0204 


0.8733 


0.7825 





Pi(u — U.42d1) 


1 X 10^ ' 


3 X 10^' 


7 X 10" ' 


1 X lO"'' 


2 X 10"" 




P2iu = 0.4261) 


2 X 10"** 


1 X 10^' 


3 X 10"' 


5 X 10" ' 


2 X 10"" 




AE{u = 0.6818) 


1.6249 


1.3121 


1.1128 


0.9907 


0.9134 


0.8463 


pi(u = 0.6818) 


9 X lO'** 


3 X 10" ' 


5 X 10"' 


1 X 10"'' 


2 X 10"" 


3 X 10"" 


P2{u = 0.6818) 


2 X lO"** 


8 X 10"** 


1 X 10"' 


2 X 10" ' 


2 X 10"' 


7 X 10"' 



TABLE I: Energy width AE, truncation errors pi (first DMRG step), p2 (second DMRG step) for u = 0, u = 0.4261 (near the 
quantum critical point), and for u — 0.6818 in the Hubbard lattice when mi = 512 and m2 — 96 states are retained. 
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FIG. 2: Scaled spin gap in the quasi-lD Heisenberg model 
as function of Jy. 



The accurate location of the critical point is done by 
plotting the product L~^^ as function of the interaction 
driving the transition. ^ is the correlation length. This is 
because at the transition, L~^^ is independent of L^. For 
the gap the function L~^S^ translates to L~^A"^, where z 
is the dynamical exponent. Near the the quantum critical 
point, the product L^A is given by a universal function, 

L^A = f{{g - gc)Ll^n, (7) 

where g is a generic coupling driving the transition, gc is 
its magnitude at the quantum critical point, and v is the 
correlation length critical exponent. 



In the second step, we targeted ritarg — 3 charge 
sectors Ne,Ne ± 1 with the lowest Sz- The reduced 
density-matrix was formed by attributing an equal weight 
Wfe = 1/3 for each of fc = l,ntarg states. We kept 
ms2 — 96 states such that the width of the retained 
states, AE » ty,td for td = ty = O.Obtx. AE is dis- 
played in Table HI For these parameters, the truncation 
error during the second step was such that p2 ^ pi for 
systems of up to L^x Ly = 32x 33 when three superblock 
states were targeted. We empirically chose ms2 such that 
AE/ty = 10. For this ratio, we can accurately reproduce 
the exact result at u = 0. 



IV. FINITE-SIZE SCALING 
A. General concepts 

We apply finite-size scaling^- to analyze the results on 
the charge gap A. The procedure is simple. We accu- 
rately compute A in order to locate the quantum critical 
point. We then collapse the data using the exponents i' of 
known universality classes in order to find the class cor- 
responding to the Mott transition. We emphasize that in 
this procedure there is no extrapolation or external pa- 
rameter besides the data and the exponent of the chosen 
universality class. 



B. Application to coupled Heisenberg chains with 

S = 1 

In Fig[^ we illustrate the finite-size analysis that we 
apply below to weakly coupled Heisenberg chains with 
S = 1. The model which was studied in Refi^ is given 
by the Hamiltonian, 

Hs = Jx ^ia:,iy^ia: + l,iy + "^V ^ ,iy^ ,iy + l ■ (8) 

In the model ((SI), there is transition from a magnet- 
ically disordered ground state, the Haldane gap phase, 
to a magnetically ordered ground state which is induced 
by the transverse coupling Jy. This transition has been 
studied by the quantum Monte Carlo method^^. In 
this transition z = 1, and it belongs to the universal- 
ity class of the 3D classical Heisenberg model, for which 
V = 0.7048^^. In FigHwe plot La^A^ as function of Jy, 
where As is the spin gap. We studied systems ranging 
from L.j; X Ly = 12 X 13 to 24 x 25. We applied peri- 
odic boundary conditions along the x-direction and open 
boundary conditions along the y-direction. At the quan- 
tum critical point Jy = Jy, L^Ag is independent of L^. 
There are small size effects for smaller systems. We thus 
included only systems larger than 16 x 17. All the curves 
L^As cross at Jy. The critical point Jy = 0.04368 was 
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FIG. 4: Error in the ground-state energy for quasi-one- 
dimensional systems as function of the hnear dimension Lx 
of the lattice. Single-particle two-step DMRG gaps versus 
exact gaps as function of La; 



fictitious universality class with ly — l.b. As expected 
from Monte Carlo simulations, the best data collapse was 
obtained for v x 0.7048 which is predicted Monte Carlo 
valued for the classical 3D Heisenberg universality class. 
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FIG. 3: A X I/a; as function of it ((a) and (b), as function of 
{Jy ~ Jy)L]J'^ for different Lx x Ly and for different univer- 
sality classes: mean-field {v = 0.5), classical 3D Heisenberg 
{v = 0.7048), 2D Ising {v = 1.0), fictitious class {v = 1.5). 



located graphically. It is in perfect agreement with the 
quantum Monte Carlo value = 0.043648(8). 

The determination of the universality class is done by 
plotting La;As as function of {Jy — Jy)Ll^'' . In Figl^l 
Lj-Ag is displayed for different values of v corresponding 
to mean-field, classical 3D Heisenberg, 2D Ising, and a 



V. RESULTS AND DISCUSSION 

A. Correlation length exponent at the Mott 
transition 

We can now confidently apply the same method to the 
Hubbard model. It has roughly the same level of diffi- 
culty as the coupled Heisenberg chain problem. First, we 
compared the two-step DMRG results with the exact en- 
ergies at u = 0. We emphasize that this test is non-trivial 
for a real-space technique such as the DMRG because in 
real space, the hopping term is non-diagonal. In Fig|4l[a), 
we show the error 6E in the ground-state energies per site 
for systems ranging from Lj; X Lj, = 12 x 13 to 32x33. The 
two-step DMRG is in very good agreement with the exact 
result; 6E < 10^® and increases relatively slowly with 
for systems x Ly < 28 x 29 and starts to grow sharply 
beyond this size. In FigHJJb), we compare the single- 
particle gap, A = ^[EoiN + l) + Eo{N - 1) - 2EoiN)], 
obtained with the two-step DMRG to the exact gap. 
The largest error for the gap was about 5 x 10^** in the 
32 x 33 systems. Since for this size the exact gap is only 
A — 0.00103, we excluded the 32 x 33 systems from the 
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FIG. 5: Quasi-particle gaps as function of Lx for two char- 
acteristic values of the interaction: (a) u = 0.2273, (b) 
u = 0.6818 for ID (circles) and quasi-lD (squares) systems. 



data used to extract the critical exponent. For the largest 
systems kept for the analysis 28 x 29, the two-step DMRG 
gap is A = 0.00895 which is to be compared to the exact 
gap A = 0.00883. The relatively large loss of accuracy 
in the gap for 32 x 33 systems follows from the sharp 
increase in 5E. 

When u ^ 0, the two-step DMRG retains the same 
level of accuracy as at w = 0. This is because, when 
the same number of states m2 is kept, the truncation 
error p remains close to that of w = as seen in Ta- 
ble HI A£' slightly increases with m, hence, the condition 
/S.E ^ ty^td is also fulfilled. Unlike the pure ID model, 
the metallic phase is expected to have a finite width in 
the quasi-lD model. In Fig. [S]we show the gap as func- 
tion of Lx for two characteristic values of the interaction 
at u = 0.2273 and u = 0.6818 for the ID and quasi-lD 
systems. There appears to be two regimes. In FiglSJa), 
for u = 0.2273 the quasi-lD gap shows a sharp decay in 
contrast to the ID gap which decays more slowly. This 
is consistent with the finite value of the ID gap and the 
presumably zero value of the quasi- ID gap in the ther- 
modynamic limit. In FiglS^b), for u — 0.6818 both gaps 
remain very close and have a finite value in the thermo- 
dynamic limit. This behavior suggests that there would 
be a quantum critical point at 0.2273 ^ Uc ^ 0.6818. We 
would like to emphasize that in Ref.^^, in ID in agree- 
ment with the exact result the DMRG yielded Uc = 0. 

We analyze our results using the language of second 
order transitions. This is justified because we did not 
see any sharp change in our data for the ground-state 



energy or the gap. Generally, in a first order transition it 
would usually be expected that the ground-state energy 
would be non-differentiable and the gap would show a 
discontinuity at the transition point. These were not seen 
in our data. The absence of a discontinuity is seen for 
instance in the behavior of Lj,A in FigjS] This justifies 
the assumption that the transition is of second order. 

As for the Heisenberg model above. In RefJ^ii^, the 
value z = 1 was predicted for the interaction induced 
Mott transition. But in the density induced transition 
the dynamical exponent is z = 2. In order to find the 
value of z, we plotted both L^jA and i^A. However, 
the rough estimate of the critical value found for L^A, 
Uc ~ 0.1705 was very inconsistent with the direct ex- 
trapolation of the data. For instance, at u = 0.2273, A 
extrapolates to 0. This allows us to rule out z = 2 as well 
as higher values of z since they yield even smaller Uc- 

We show for z = 1, La; A as function of u in FiglH] 
A first sweep of the interaction range < u < 0.6818 
in Figita) indicates that 0.4 < < 0.5. In Figitb), 
to precisely locate Uc, we concentrate in the interaction 
range 0.420 < m < 0.432, a graphical estimate yields 
Uc — 0.4255. The range of values of m for the critical anal- 
ysis Su = 0.02656mc is comparable to that used in Refi^ 
\S{J/U)\ = 0.01526(J/;7)c for the Bose Hubbard model, 
and in Ref}^ \S{t/g)\ = 0.01339(i/5)c for the Jaynes- 
Cummings-Hubbard model. (J/U) and t/g are the ratio 
of the hopping parameter over the interaction. 

As for the Heisenberg model above, we determine the 
universality class of the Hubbard model by plotting L^jA 
as function of {u — Uc)L^^'^ ■ In Fig[71 We tried different 
values of corresponding to the mean-field h' = 0.5, 3D 
XY, 2D Ising = 1.0, and a fictitious = 1.5 cases. 
For the 3D XY model, Monte Carlo values v are found 
between i/ 0.662(7) and = 0.6723^^, and with the 
bosonic Hubbard modelH and the Jaynes-Cummings- 
Hubbard modeli^ for which v = 0.6715. The experiments 
on He films are believed to yield the best estimate of v 
for the 3D XY models. Experiments have smaller errors 
than Monte Carlo simulations. For instance v was found 
to be = 0.6708(4) in Ref A = 0.6705(6) in Ref-^, 
and 1/ = 0.67095(13) in Ref,^. We used this last value to 
collapse the data for the test of the 3D XY universality 
class. 

FigUl clearly shows that the best fit to the data is ob- 
tained for 1/ = 1.0. This implies that the Mott transition 
in the Hubbard model belongs to the universality class 
of the 2D Ising model as predicted by the approximate 
mapping of Rcf.^^^. 

The 3D XY universality class for the Mott transition 
in 2D was conjectured in approximate slave-rotor anal- 
yses of the fermionic Hubbard model in Ref^^ii^. This 
work shows that the neglect of the gauge field during the 
factorization of the fermionic operators into a spinless 
boson and a charge-less spin is not justified. It should be 
noted that the 3D Ising and 3D Heisenberg universality 
class for which i/ is close to that of the 3D XY class, re- 
spectively 1/ = 0.6298(5)?£, v = 0.7048S were also ruled 
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FIG. 6: A function of u for the Hubbard model: (a) 

extendend range of u, (b) for u in the vicinity of the quantum 
critical point. 
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B. Order parameter for the Mott transition 
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The identification of the universality class of the Mott 
transition suggests the following analogy with the Ising 
transition. The weak u limit should correspond to the 
high temperature phase in the Ising model. At u = 0, 
the four possible local states, |0), | t), I i), and | ti) are 
equally probable respectively with 



no = = ni = nd 



1 

4' 



(9) 



in n-j- and hi only purely singly occupied sites are 
counted, 



(10) 
(11) 



In the opposite limit u ^ oo which corresponds to the 
low temperature phase, holes and doubly occupied sites 
are not allowed, 



no = rid = 0, 
1 



(12) 
(13) 
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FIG. 7: A X La; as function of it ((a) and (b), as function of 
{u — Uc)Ly'^ for different Lx x Ly for v corresponding to dif- 
ferent universality classes: mean-field {v — 0.5), 3D classical 
XY {u = 0.67095), 2D Ising {u = 1.0), fictitious {u = 1.5). 



the local possible states have shrinked from 4 to 2 due 
to the Z2 Ising symmetry breaking. This is in con- 
trast to the slave-rotor analyses where the Fermi liq- 
uid is regarded as the ordered phase. The isomorphism 
SU{2)/Z2 = 5*0(3) implies that in principle after the 
Mott transition, the effective spin Hamiltonian, obtained 
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by projecting out the empty and doubly occupied states, 
should retain the full spin rotational symmetry. The 
eventual spin long-range order will depend on the cou- 
plings present in the effective Hamiltonian. 

The natural order parameter M for the Mott transition 
should thus be given by the average number of singly 
occupied sites minus the number of doubly occupied and 
empty sites, 

M = {n^ + n^-nd-riQ). (14) 

Thus, 

M={n- Ana). (15) 

For u = 0, M = and for u = oo, M = 1. It should 
be expected that for u < Uc, M — 0. But this is not true 
for finite systems. Because of the finite size gap, finite 
systems are always 'ordered', thus M will always have a 
finite value for a finite system even when u < Uc- Since 
Ud is a local quantity, it changes very slowly with system 
sizes. This means that very large systems are necessary 
to extrapolate accurately to its thermodynamic value. 
In the two-step DMRG approach , it is more judicious to 
calculate the correlator. 



for the middle chain. M is shown in FiglS] for a 24 x 25 
system. The curve of M has the usual form of an order 
parameter curve. However in the vicinity of the quantum 
critical point, because of the use of open boundary condi- 
tions, the data are strongly affected by the 2D remnant 
of Friedel oscillations. Convergence is very slow even 
with this definition of M. It can be seen that the value 
of M is still appreciable at the quantum critical point 
u = 0.4255. Significantly more work will be necessary in 
order to reliably extract the order parameter exponent 



more realistic studies of the detailed and well controlled 
analysis of the Mott transition. The critical behavior of 
the quasi-lD model Hubbard model is found to belong 
to the universality class of the 2D Ising model. The fact 
that the transitions in the quasi-lD Heisenberg and Hub- 
bard models belong to the universality classes of their 
isotropic counterparts shows that despite the restriction 
of the two-step DMRG method to highly anisotropic 2D 
models, it is nevertheless very useful for the understand- 
ing of the physics of isotropic 2D systems. 

We did not discuss the spin degrees of freedom. They 
are expected to be gap-less in either side of the Mott 
transition. In the insulating phase, in the strong coupling 
limit U 3> tx,ty,td, the anisotropic frustrated Hubbard 
model is equivalent to the anisotropic Ji — J2 model with 
Jx = t^/U , Jy = ty/U, and Jd = td/U, where Jx, Jy, and 
Jd are respectively the exchange parameter in the longi- 
tudinal, transverse, and diagonal directions. Our choice 
ty = td implies that the ground state will be magnetically 
ordered with the momentum q = (7r,0). For intermedi- 
ate U, in the Mott insulator phase, double occupation is 
not negligible straighforward mapping to the Heisenberg 
model is not valid. However, the charge gap opening im- 
plies spin-charge separation. Thus even in this case, the 
effective low-energy Hamiltonian should be Heisenberg 
like, albeit with non-trivial exchange parameters. Mag- 
netic long-range order should be expected. However, a 
gap-less spin-liquid ground state with a spinon Fermi sur- 
face as suggested in Ref.^t?^ jg ^igQ possible. 
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VI. CONCLUSION 

In this paper, we used the two-step DMRG to analyze 
the finite size behavior of the quasi-particle gap in the 
ground-state Mott transition in the quasi-lD Hubbard 
model. We chose a non-bipartite lattice to avoid the issue 
related to the possible nesting induced Slater transition. 
We studied systems ranging from 12 x 13 to 32 x 33. 
We were able to find the universality class of the Mott 
transition in an un-biased calculation. 

In contrast to the pure ID model, we find that the 
quasi-lD models displays a genuine Mott transition at 
a finite critical interaction. Moreover, the quasi-lD so- 
lution does not have the pathologies of the infinite di- 
mensional solution. It could thus serve as a basis for 




u 



FIG. 8: Order parameter M of the Mott transition as function 
of u for a 24 X 25 system. 
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